#include <UnitTest++.h>
#include <Utility/CheckPoint.hpp>
#include <3rdParty/TaucsWrapper.hpp>
#include <zMat.hpp>

using namespace zzz;
#ifdef ZZZ_LIB_TAUCS
SUITE(TaucsTest)
{
  TEST(SymmetricTest)
  {
    // 1    0.5   0     0
    // 0.5  1     0.5   0
    // 0    0.5   1     0.5
    // 0    0     0.5   1
    zMatrix<double> a(Zerosd(4,4));
    a(0,0)=1;
    a(0,1)=0.5;
    a(1,0)=0.5;
    a(1,1)=1;
    a(1,2)=0.5;
    a(2,1)=0.5;
    a(2,2)=1;
    a(2,3)=0.5;
    a(3,2)=0.5;
    a(3,3)=1;
    zVector<double> b(4);
    b(0)=1;
    b(1)=2;
    b(2)=3;
    b(3)=4;

    zVector<double> x(Inv(a)*b);

    // TaucsSolver accepts and only accepts symmetric matrix with only lower part filled
    zSparseMatrix<double> sa(Tril(a));
    TaucsCholeskySolver solver;
    CHECK(solver.Creater(sa,"none"));
    zVector<double> sx(4);
    CHECK(solver.Solve(sx.Data(), b.Data()));

    for (zuint i = 0; i < x.Size(1); i++) {
      CHECK_CLOSE(x[i], sx[i], EPSILON);
    }
  }

  TEST(NonSymmetricTest)
  {
    //  0 2 1
    //  2 0 1
    //  0 1 2
    //  3 0 0
    zMatrix<double> a(Zerosd(4,3));
    a(0,1)=2;
    a(0,2)=1;
    a(1,0)=2;
    a(1,2)=1;
    a(2,1)=1;
    a(2,2)=2;
    a(3,0)=3;

    zVector<double> b(4);
    b(0)=2;
    b(0)=3;
    b(0)=4;
    b(0)=5;

    zVector<double> x(Inv(Trans(a)*a)*Trans(a)*b);

    // TaucsSolver accepts and only accepts symmetric matrix with only lower part filled
    zSparseMatrix<double> sa(Tril(Trans(a)*a));
    TaucsCholeskySolver solver;
    CHECK(solver.Creater(sa,"none"));
    zVector<double> atb(Trans(a)*b);
    zVector<double> sx(4);
    CHECK(solver.Solve(sx.Data(), atb.Data()));

    for (zuint i = 0; i < x.Size(1); i++) {
      CHECK_CLOSE(x[i], sx[i], EPSILON);
    }

  }
}
#endif // ZZZ_LIB_TAUCS